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Abstract. We study the formation of (quasi-)coherent matter waves emerging from a 
Mott insulator for strongly interacting bosons on a one-dimensional lattice. It has been 
shown previously that a quasi-condcnsate emerges at momentum fccond = 7r/2a, where 
a is the lattice constant, in the limit of infinitely strong repulsion (hard-core bosons). 
Here we show that this phenomenon persists for all values of the repulsive interaction 
that lead to a Mott insulator at a commensurate filling. The non-equilibrium dynamics 
of hard-core bosons is treated exactly by means of a Jordan- Wigner transformation, 
and the generic case is studied using a time-dependent density matrix renormalization 
group technique. Different methods for controlling the emerging matter wave are 
discussed. 
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1. Introduction 

Cold atoms in optical lattices currently constitute one of the most flexible and, 
hence, most promising ways to investigate areas of many-particle physics for which no 
established knowledge prevails. One such area is that of strongly correlated systems, for 
which neither the ground-state nor excited states can be easily described and for which 
no generic analytic tools with predictive power are presently available. The situation is 
even more difficult for non-equilibrium processes, where, in general, a full knowledge of 
the eigenstates of the system is needed in order to properly describe its dynamics out 
of equilibrium. 

Recent experimental work has succeeded in producing one of the hallmarks of strong 
correlations, namely, a Mott insulator 0I2]- Remarkably, this was achieved with bosons, 
a situation not easily encountered in condensed matter physics. Even the limit of very 
strong interactions, where the bosons can be considered to be impenetrable particles 
(hard-core bosons or a Tonks- Girardeau gas could be reached experimentally 
in |3] . Such a limit is attainable in elongated traps for large positive three-dimensional 
scattering lengths, at low densities, or with very strong transversal confinement [HI 13 IE] ■ 

On the other hand, the study of the non-equilibrium dynamics of quantum gases was 
instrumental for understanding their properties [H]. More recently, a controlled study 
of the evolution of one-dimensional Bose gases prepared initially out of equilibrium was 
performed ^Uj, opening up the possibility of investigating the role of integrability in the 
relaxation dynamics of a many-body system In the case of hard-core bosons in one 
dimension, the evolution of a system with an arbitrary initial state can be described 
theoretically in an exact way ^21 ^1 exact mapping onto free fermions, the 

Jordan- Wigner transformation 

We concentrate here on the out-of-equilibrium evolution of a one- dimensional Bose 
gas with strong interactions whose initial state is a Mott insulator. It was previously 
shown that, in the hard-core limit, an initial Fock state develops quasi-long-range 
correlations when the bosons are allowed to evolve freely on a lattice ^21 E]- In 
particular, the momentum distribution function = (rik) develops sharp peaks at 
momenta k = ±7r/2a, where a is the lattice constant. An examination of the one- 
particle density matrix shows that, after the formation of the peaks in rik, it decays 
as \l\fx at long distances, as it does for hard-core bosons in equilibrium [T7j . 
demonstrating that, in fact, the peaks in signal the emergence of quasi- coherence 
at a finite wavevector. A detailed picture of the (quasi-) coherent part is obtained 
by examining the lowest natural orbital (NO), i.e., the eigenvector of the one-particle 
density matrix corresponding to the largest eigenvalue. Once the peaks in form, 
the NO evolve at a constant velocity ujvo = ±2at//i, where t is the nearest neighbor 
hopping amplitude, without appreciable change in their form. These are the maximal 
group velocities on a lattice with a dispersion = —It cos ka. The process of formation 
of the quasi-condensate is also characterized by a power law. The population of the 
quasi-condensate increases in a universal way as ~ 1.38^tr//l, as a function of the 
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evolution time r, independently of the initial number of particles in the Fock state. The 
time Tm at which the maximal occupation of the NO is reached depends linearly on the 
number of particles Ni, in the initial Fock state, and is given by Tm = 0.32Ni,h/t. 

The appearance of quasi-condensates at = ±11 /2a can be understood on the 
basis of total energy conservation. Given the dispersion relation of hard-core bosons 
on a lattice, since the initial Fock state has a flat momentum distribution function, 
its total energy is Ej- = 0. If all the particles were to condense into one state, it 
would be to the one with an energy Ik = Et/N. Taking into account the dispersion 
relation = —2t cos ka, Ik = corresponds to k = ±7r/2a. Actually, since there is 
only quasi-condensation in the one-dimensional case, the argument above applies only 
in that the occupation of a given state is maximized. In addition, the minimum in 
the density of states at these quasi-momenta strengthens the quasi-condensation into a 
single momentum state. 

Since hard-core bosons can be treated exactly [U [121 ESI Ell 113 , they are extremely 
well-suited for a theoretical study of nonequilibrium dynamics because large systems 
(with hundreds to thousands of bosons) can be examined over long times. However, 
the experimental investigation is hampered by the quite stringent requirements for the 
realization of such systems. We therefore consider here the case of finite interactions, 
modeled by the one-dimensional Hubbard model 

^ = E (^l^'+l + ^•^•) +JJ2''^{n^-l), (1) 

i i 

where h\ and are bosonic creation and annihilation operators, respectively, and 
rii = h\h^ is the density operator. The hard-core limit corresponds io U 00. The 
value at which the Mott insulator appears has been estimated as Uc/t ~ 3.5 in one 
dimension for a commensurate density n = 1 |lH|. Hence, all the cases considered here 
correspond to U > Uc- 

As in the hard-core case, we start with bosons in a Mott-insulating state spread over 
several lattice sites and monitor the free expansion on a lattice. For the time evolution 
of the system, it is, in principle, necessary to treat the whole Hilbert space of the system, 
restricting exact treatments to extremely small systems. Instead of fully diagonalizing 
the Hamiltonian matrix, efficient iterative eigensolvers such as the Lanczos or the Jacobi- 
Davidson procedure are commonly used ^011201, enabling one to treat somewhat larger 
systems. 

Unfortunately, the Lanczos method is limited by the exponential growth of the 
Hilbert space as a function of the number of degrees of freedom. Therefore, a 
more efficient way of representing the relevant subspaces for the time evolution is 
needed. Recent progress in this direction was achieved by extending the density matrix 
renormalization group (DMRG) method pT| l^ to treat the time evolution of correlated 
systems [231121112211211123123, leading to the so-called t-DMRG. 

Here we apply the t-DMRG to study the expansion of soft-core bosons out of a 
Mott insulator and compare it to the hard-core case. It will be shown that the essential 
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features obtained with hard-core bosons are preserved, while new control possibilities 
are opened by tuning the strength of the interaction U . We will also show that 
further control of the momentum of the emerging quasi-condensates can be achieved 
by introducing a superlattice, which is obtained by superimposing an extra periodic 
potential onto an already existing lattice potential. Such systems have been recently 
reahzed experimentally with ultracold gases trapped on optical lattices jSOJ EI], and 
have been studied theoretically by various mean field approaches [32], quantum- Monte 
Carlo simulations and exact diagonalization pM] , 

In Sec. 121 we discuss the theoretical treatment of the time evolution of a general 
quantum system both within the framework of the Lanczos method and of the t-DMRG. 
The results are shown in Sec. OJ where the evolution for parameters in the range 
6 < f//t < 40 are considered. As in the hard-core case, the momentum distribution 
function Uk displays maxima at finite wavevectors. However, those wavevectors are 
displaced to lower values of k as the strength of the interaction is reduced. In Sec. El 
we study how the introduction of a superlattice allows further control of the emerging 
quasi-condensates. For these systems we restrict our analysis to the hard-core regime. 
Finally, a concluding discussion is given in Sec. El 

2. Time evolution of many-body quantum systems 

The general solution of the Schrodinger equation for a many-body system can only be 
given in a formal way, that for a time-independent Hamiltonian is 

|^(r)) = e-^^-|^(r = 0)), (2) 

where H is the Hamiltonian determining the evolution over a time r from some initial 
time r = 0. For sufficiently small systems, full diagonalization of H is possible, and a 
knowledge of all the eigenvalues and eigenstates allows for an exact determination of the 
evolved state. However, for a general many-body system, the size of the Hilbert space 
grows exponentially with the number of degrees of freedom, restricting the system size 
essentially to tens of atoms. More efficient ways of determining the time evolution are 
discussed in the following subsections. 

2.1. Lanczos method 

Here we focus on the Lanczos procedure, which can be generalized in a straightforward 
way so that the time evolution of the system can be computed without calculating all 
eigenstates. 

The basic idea is to expand the time-evolution operator, 

I ^(r + Ar) ) = e-^^^- | ^(r) ) ^ I ^{r) ) , (3) 

n=0 

and to focus on the set of states {| ?/'(r) ), | ^/'(r )),... , iJ"* I ^(t) )}, which spans the 
so-called Krylov subspace. 
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The key idea of the Lanczos method is to obtain a basis by orthogonahzing the 
vectors of the Krylov subspace only with respect to the previous two elements of the 
set, leading to the recursion relation 

I Vi+i ) = H\ Vi ) -ai\ Vi ) - I Vi_i ) (4) 

{Vi\H\Vi) 2 {Vi\Vi) 

with ai = — — ■ — — , (3i = ■ (5) 

{Vi\Vi} {Vi^i\Vi_i} 

for the vectors \vi) of the Lanczos basis. The projection of the Hamiltonian onto this 
basis set leads to a tridiagonal matrix = V^HVm, where is a rectangular matrix 
containing the Lanczos vectors as column vectors. In this way, the Hamiltonian can 
be efficiently diagonalized. 

Using this approach, an approximation for the time evolution of a given state 
\iP{t)) at time r over a small time interval At can be given. For a time- independent 
Hamiltonian it reads 

I ^(r + At) )„ = V^^-^^-^^Vl I ^{t) ) . (6) 

Remarkably, an exact bound can be given for the error in the approximation 



^{t^At) )- I V'(r+Ar) )approx II < 12 exp 



16m 



( epA 



T 



\ Am 



(7) 



valid for m > pAr/2, where p is the width of the spectrum of H and m is the dimension 
of the Krylov subspace. For pAr ^ 1 and m sufficiently large (of the order of 10), the 
formula above shows almost exponential convergence. 



2.2. Adaptive time evolution with the DMRG 

The basic idea of the density-matrix renormalization group method is to represent one 
or more pure states of a finite system approximately by dividing the system in two and 
retaining only the m most highly weighted eigenstates of the reduced density matrix of 
the partial system. In combination with the numerical renormalization group approach 
(NRG) developed by Wilson jHEj and the superblock algorithms developed by White 
and Noack j2ZI, this leads to a very powerful and efficient tool for the investigation of 
one-dimensional strongly correlated quantum systems on a lattice. Here we only give 
a rough sketch of the method and refer to recent reviews [201 iSHl 121] for a detailed 
description. 

As depicted in Fig.d the key steps are to increase the number of degrees of freedom 
of the partial system by adding sites, then to decrease the number of degrees freedom 
by retaining states below a cutoff. In this way, the method carries out a renormalization 
group procedure closely related to Wilson's NRG. 

In the first step of the algorithm, a site is added to one of the subsystems, with 
its Hamiltonian exactly represented. The Hamiltonian of the subsystem is usually 
represented in an efficient reduced basis built up from the m most important eigenstates 
of its reduced density matrix. Note that the basis is incomplete due to the truncation. A 
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RG-step 1 : Increase number of degrees of freedom: 
Add exact site to old system block 



Obtain \W> 



Obtain p 
Diagonalize p 



New basis: eigenstates of p 
RG-step 2: Decrease number of degrees of freedom: 
Cutoff after m states 



Transform system block 
into new basis with only m states 



Figure 1. Sketch of the lattice and flowchart of the DMRG iteration scheme. 
The lattice is shown in the usual "superblock" configuration, where the left part of 
the lattice is the subsystem which is used to compute the basis of density matrix 
eigenstates. At the 'dividing' bond, two "exact" sites are added; the "sweep" proceeds 
from left to right. The flowchart at the right shows the relevant steps of the DMRG 
procedure as described in the text. 



measure of the error e introduced by the cutoff is the discarded weight, which measures 
the total weight of the discarded states 

m 

where Xj is the j*^ eigenvalue of the reduced density matrix. After convergence is 
reached, one typically obtains values e < 10~^ with m < 1000. Although one is only 
retaining a tiny fraction of the total Hilbert space of the system (which already for small 
systems can reach a dimension of several million), the desired observables can thus be 
obtained with high accuracy. 

In the second step of the iteration, the states one is interested in are obtained. 
These states are called "target states". In the original ground-state algorithm, these 
are the ground state and the few lowest lying excited states of the system, which are 
obtained by diagonalizing the Hamiltonian of the total system, e.g., by carrying out the 
Lanczos diagonalization algorithm described in the introduction. However, states other 
than the eigenstates of the Hamiltonian may be obtained in this step. This flexibility 
is crucial for the time-evolution algorithms, in which the time-evolved state is obtained 
by other means than solving an eigenvalue problem. 

In the third step, the new effective basis is obtained by diagonalizing the reduced 
density matrix of the extended subsystem, given by 

Psubsystcm = Truest rii \ tpi) { tpi \^ , ^ rij = 1 , (9) 

where the sum goes over all target states. In step four, only the m eigenstates with 
the largest eigenvalues are kept. The operators needed to represent the Hamiltonian 
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of the subsystem, to form the pieces of the Hamiltonian connecting subsystems, and 
to calculate observables are transformed into this new reduced basis. This effective 
Hamiltonian of the subsystem is now the starting point for step one of the next iteration. 
In this way, every step improves the accuracy of the obtained eigenstates and energies 
by improving the reduced basis used for the representation of the target states. 

DMRG iteration schemes are usually divided into two classes. In the so-called 
infinite- system algorithm, the system grows at each step. This can be used to build 
up the system up to a desired lattice size. In the finite- system algorithm, the size of 
the lattice is fixed, and the "dividing bond", i.e., the position at which the system is 
cut in two parts, is moved from the right end of the lattice to the left end and back 
(other variations are possible). This is called a "sweep". In order to obtain the ground 
state (and optionally the lowest lying excited states) with a high accuracy, the sweeps 
are iterated until convergence is reached. The calculation can be significantly sped up 
if the diagonalization in step 2 of the DMRG procedure is started with a good initial 
guess for the wave function. Such an initial guess can be constructed using the so-called 
"wave function transformation", which approximately transforms the wave function 
obtained from the previous finite-system step into the basis of the current superblock 
configuration. As we will see next, the wave function transformation also plays a key 
role in the adaptive t-DMRG schemes. 

The main difficulty in calculating the time evolution using the DMRG is that the 
restricted basis determined at the beginning of the time evolution is not able, in general, 
to represent the state well at later times [21] because it covers a subspace of the total 
Hilbert space which is not appropriate to properly represent the state at the next time 
step. Since both the Hamiltonian and the wave function | iI){t) ) at time r are represented 
in an incomplete basis, the result for the next time step | ?/;(r + Ar) ) will have additional 
errors because the reduced basis is not an optimum representation for this state. In 
order to minimize these errors, it is necessary to form a density matrix whose m most 
important eigenvectors are "optimal" for the representation of the state | iP{t) ), as well 
as for I ipij + Ar) ) in the reduced Hilbert space. The most straightforward approach 
is to mix all time steps | V'(Ti) ) into the density matrix [211 122] • However, this can be 
extremely costly computationally. A more efficient way is to adapt the density matrix 
at each time step. 

An approach for adaptive time evolution based on the Trotter-Suzuki |23] 
decomposition of the time-evolution operator was developed in Refs. [201 123 I2H1 • The 
idea is to split up the time-evolution operator in local time-evolution operators Ui acting 
only on the bond I. For lattice Hamiltonians containing only terms connecting nearest- 
neighbor sites, this is easily obtained using the Trotter-Suzuki decomposition, which in 
second order is given by 

^-iArH ^ g-jArHcvGn/S ^-iArHodd g-jArJ^cven/S (10) 

Here -ffcvcn and i^odd is the part of the Hamiltonian containing terms on even and odd 
bonds, respectively. Since each bond term Hi within i/cvcn oi H^dd commutes, e"*^"^^ 
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Apply local U/ at the bond dividing 
the system and the environment block 










Obtain m density-matrix basis states 








Shift the 'dividing bond' by one lattice site, 
using the 'wave function transformation' 



Add the following states to the density matrix: 
IV(t:)> IV(T+Ax/ii)>l¥(x+2AT/ii)>- |v|/(x+At)> 



Obtain m density-matrix basis states 



Shift the 'dividing bond' by one lattice site 
using the 'wave function transformation' 



Figure 2. The flowcharts of the t-DMRG schemes described in the text. On the left, 
the flowchart for the Trotter variant is sketched. On the right, the scheme used for the 
Lanczos variant is shown. 



can then be factorized into terms acting on individual bonds. As depicted in Fig. ^ in 
the DMRG procedure usually two sites are treated exactly, i.e., the entire Hilbert space 
of the two sites is included. The Trotter variant of the t-DMRG exploits this feature by 
applying Ui = e"*^"^^' at the bond given by the two "exact" sites. In this way, the time- 
evolution operator has no further approximations other than the error introduced by 
the Trotter decomposition. In particular, the error introduced by the cutoff is avoided. 
The wave function of the lattice is then updated by performing one complete sweep over 
the lattice and applying Ui at the "dividing bond" . In this way, only one wave function 
must be retained and it is possible to work with the density matrix for a pure state. 
The flowchart is sketched in Fig. |2l 

However, the method is restricted to systems with local or nearest-neighbor terms 
in the Hamiltonian. A more general basis adaption scheme aims at adapting the density 
matrix basis by approximating the density matrix for a time interval jlO] , 

T+Ar 

PAr= J |^(r'))(^(r')|rfr'. 

T 

The integral is approximated by adding a few intermediate time steps within the time 
interval [r, r + Ar]. In Ref. jlOj, the intermediate time steps are obtained using a Runge- 
Kutta integration scheme and using 4 or 10 intermediate time steps. Here we instead 
obtain the intermediate time steps using the Lanczos method described in Sec. 12.11 This 
can be done easily because the Hamiltonian of the system is usually constructed anyway 
in the DMRG scheme. Within the restricted basis, the Lanczos iteration, Eq. can 
then be performed, leading to the desired intermediate time steps computed using Eq. 
(jH)). Using this approach, we find that, if the time step is small enough, it is sufficient 
to retain only the target state | ip{T + Ar) ), so that one can work with a pure-state 
density matrix like in the Trotter approach. For larger time steps, it is important to mix 
at least the states \iP{t)) and \ijj{T + Ar) ) into the density matrix. We find that it is 
sufficient to perform only one half-sweep in order to adapt the restricted basis, which is 
the minimum requirement for updating the basis on the complete lattice. The flowchart 
of this approach is sketched in Fig. [21 With this approach, it is, in principle, possible 
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Figure 3. Error analysis of the t-DMRG obtained by comparing the results for the 
density (n^ ) and of the momentum distribution Uk of an initial Fock state with 10 
hard-core bosons on a lattice with 50 sites with exact results obtained using the Jordan- 
Wigner transformation. The t-DMRG results were calculated using the Trotter variant 
(second order) of the method, keeping up to m — 200 basis states. 

to treat more general Hamiltonians, as long as they can be treated accurately using 
the DMRG. However, due to the fact that, in general, one cannot work with a pure 
state density matrix, the dimension m of the restricted basis needed to obtain a certain 
discarded weight during the time evolution is larger compared to the Trotter approach 
described above, making this variant slower. The errors in both adaptive schemes for 
intermediate and long times are comparable. For the problem at hand, we therefore 
choose the Trotter variant (in second order) of the adaptive time-dependent DMRG. 

In this work, we control the error during the time evolution by fixing the discarded 
weight and varying the number of basis states kept. By keeping a maximum of m = 800 
density matrix eigenstates, we obtain discarded weights smaller than 5 ■ 10"^ during the 
time evolution. 

The initial state has zero density over a wide region in the system. This may lead 
to difficulties for the DMRG. In order to control this, we compare the time evolution 
of an initial Fock state of hard-core bosons obtained using the t-DMRG with the exact 
results from the Jordan- Wigner transformation. The maximum error as a function of 
time is plotted in Fig. IHl As shown in this figure, the maximum deviation from the 
exact results is smaller than 0.01 for all times considered. 

3. Free expansion of soft-core bosons from a Mott insulator 

We consider here the free expansion of interacting bosons described by the 
Hamiltonian on a one-dimensional lattice with L sites, lattice constant a and open 
boundary conditions. In the cases considered here, we take A'';, = 20 and L = 60. Due to 
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memory limitations, it is not possible to allow for all possible occupations of a given site. 
In general, the maximal number of bosons per site needed to have an accurate description 
of the system increases as U decreases. In all the cases treated here a maximum of three 
bosons per site was sufficient. Even at the smallest interaction studied here {U/t = 6), 
no appreciable difference was observed when the cutoff was changed from 3 to 4 bosons 
per site. Since the system becomes more dilute in the course of the free expansion, the 
limitation in the number of bosons per site becomes even less important at later times. 
The time sequences shown are all limited to times shorter than the time it takes the 
matter wave to reach the boundary of the system. 
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40 50 




Figure 4. Comparison between hard-core bosons (filled black symbols) and soft-core 
bosons (unfilled colored symbols) with U/t = 40 for (a) density and (b) rife at time 
r = (V), T = 2.52 (□), r = 4.98 ((j), and r = 7.5 (A) in units of 1/t. 



Figure |3] shows a comparison of the density (Fig. IH^a)) and the momentum 
distribution function (Fig. Ilfb)) for a system with U/t = 40 with hard-core bosons 
at four different times. At such high values of the interaction, it is expected that 
the particles behave as hard-core bosons. In fact, there is no noticeable difference in 
the density profiles at any time. However, the momentum distribution functions at 
r = show clear differences. While hard-core bosons are equally distributed over all 
momenta, the soft-core case has a maximum at /c = 0, showing that even in a Mott 
insulator, the fluctuations of the number of particles at each side populate that state 
preferentially. Nevertheless, after the Mott insulator is allowed to expand, the difference 
between both systems becomes barely noticeable. Already at the second time shown 
in Fig. m where a Mott plateau still exists at the center of the cloud (Fig. Ufa)), the 
momentum distribution functions of the hard- and soft-core bosons are very close to 
each other. This is expected because the constraint of a hard-core should become less 
relevant when the system is diluted. 
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Ixj - Xjl/a 




Ixj - Xjl/a 



Figure 5. (a) Modulus of the one-particle density matrix pij versus distance | xt—Xj \ 
in units of the lattice constant a at r = 0. (b) The same quantity at time r = 7.38, 
where the peaks in Uk are fully developed. Crosses (U=20 E) correspond to the one- 
particle density matrix in equilibrium for the same number of particles and system 
size. The black line corresponds to f{x) ~ I/a/I Xi — xj \. 



In order to see the establishment of coherence exphcitly, we examine the spatial 
behavior of the one-particle density matrix. We would expect a change from an 
exponential decay in the Mott-insulating state to a power law behavior if (quasi-) co- 
herence emerges. Figure |31 shows the spatial behavior of the one-particle density matrix 
both at time r = (Fig. Eta)), when bosons are in a Mott-insulating state, and at 
time r ~ 7.5 (Fig. Efb)), when the peaks around k = ±7r/2a are well established. 
Figure Et a) shows the spatial dependence of the one-particle density matrix measured 
from the center of the bosonic region in a semi- logarithmic plot for different values of 
U. As expected for a Mott insulator, an exponential decay with a correlation length 
that shortens as U is increased is observed. Figure |31^b) shows the decay of the one- 
particle density matrix on a log-log scale at a time long enough so that the peaks around 
k = ±7r/2a are fully developed. The evaluation was made in the part of the system where 
the lowest NO is appreciable, i.e. in the region of the system where a well developed 
quasi-condensate can be expected. Figure IHl shows the spatial dependence of the lowest 
NO at the same time as in Fig. Efb). The correlations in Fig. |31^b) were measured from 
the site Xj = 37a (a position where the NO is well developed) and with Xi > Xj. For 
comparison, we superimposed the one-particle density matrix for U/t = 20, A^";, = 20, 
and L = 60 in equilibrium, where, due to the lower density with respect to the initial 
state in Fig. Efa), a quasi-condensate exists. Over the distances where the NO has an 
appreciable value, no difference with the corresponding quantity in equilibrium can be 
noticed. It can be clearly seen that the one-particle density matrix has developed a 
power-law decay (the same as the one in the system in equilibrium) at the later time. 




with a power that approaches the one of hard-core bosons. Unfortunately, the expansion 
of the cloud and the total system size are not as large in the soft-core (as seen in the 
departures from the power-law due to finite-size effects) as in the hard-core case, so that 
the exponent of the power-law decay cannot be as accurately determined. Nevertheless, 
it is clear that a change from an exponential to a power-law decay takes place, indicating 
that a quasi-coherent matter wave has developed. 

Finally, we discuss the behavior of the expansion for smaller values of the interaction 
U. Although at U/t = 40 the momentum distribution function closely follows the shape 
of Uk of hard-core bosons, a tiny asymmetry can be seen in Fig. IHb) around the peaks 
ai k = IT /2a. Such an asymmetry indicates that the maximum of Uk is not exactly 
ai k = 71 /2a, but is shifted slightly. A more detailed analysis for 6 < f//t < 40 is 
presented in Fig. [7| Figure [7||^a) shows Uk around k = 7i/2a with the data points from 
t-DMRG denoted by symbols with spline interpolations between them. It is clearly seen 
that the maximum of is displaced to smaller momenta as U decreases. The spline 
interpolation allows for a better determination of the maxima in n^, since a denser set 
of fc-points corresponds to having a much longer lattice in a physical realization. On 
the other hand, the actual set of /c-points in the t-DMRG simulation corresponding to 
L = 60 is dense enough to allow for a smooth interpolation without introducing artefacts 
due to the spline procedure. Figure Ulh) shows the location of the maxima of as a 
function of U in units of 2a/7r, giving a guide for a fine tuning of the wavelength of the 
matter wave via the interaction strength. 

The results in this section show that the main feature found for the case of hard-core 
bosons ^21 1^ , namely, the emergence of a quasi- coherent matter wave from a Mott 
insulator, persists under more general conditions. The wavelength of the matter wave is 
determined primarily by the underlying lattice, on which the expansion takes place, as 
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Figure 7. (a) The momentum distribution n^. at time t = 4.5/t for different values of 
U/t. The symbols correspond to t-DMRG results on a lattice with L = 60, while the 
lines in the respective colors are spline interpolations, (b) The position of the peaks 
of the interpolated splines in nj. as a function of U/t. 



in the case of hard-core bosons. Additionally, finer tuning is possible by regulating the 
interaction strength of the bosons. This means that on an optical lattice the wavelength 
of the corresponding laser beam would essentially determine the momentum of the 
matter wave and a fine tuning can be reached by varying its intensity. The next section 
discusses further control possibilities by introducing more complex structures in the 
optical lattice. 

4. Expansion in a superlattice 

So far we have studied how finite values of the on-site interactions between bosons 
modify the behavior already known in the hard-core limit. The general finding has 
been that the physics is similar, and that an emergence of quasi-condensates can be 
obtained experimentally for a wide range of finite repulsive interactions. In this section, 
we analyze how the introduction of a superlattice potential allows a further degree of 
control over the system and can lead to a richer momentum distribution of the expanding 
cloud of bosons. We will restrict the analysis to the hard-core limit keeping in mind its 
relevance to the soft-core regime. 

In the presence of a superlattice potential, the hard-core boson Hamiltonian 
becomes 

H = -^Y1 (^l^^+i + +aY,cos (11) 

i i 
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with the additional on-site constraints 

br = b^ = 0, {6„6l} = l, (12) 

which exclude double or higher occupancy. The bosonic creation and annihilation 
operators at site i are denoted by b] and b^, respectively, and the local density operator 
by rii = b\b^. The brackets in Eq. ()12|1 apply only to on-site anticommutation relations; 
for i 7^ j, these operators commute as usual for bosons; = 0. In Eq. (fTTj). the 

hopping parameter is denoted by t and the last term represents the superlattice potential 
with strength A and I sites per unit cell. 

Using the Jordan- Wigner transformation 



i-l i-1 

b\ = flX{e-'-^'^f^ h = X{e'-^y^U, (13) 

P=l 13=1 

one can map the HCB Hamiltonian onto that of noninteracting spinless fermions, 

HF = -tY. + H.c.) + A 5^ cos ^nf , (14) 

i i 

where // and /j are the creation and annihilation operators for spinless fermions at site 
z, and n{ = flf^ is the local particle number operator. For periodic systems with 
lattice sites, one needs to consider that 

^1^^ = -/i/tv exp ^^TT^ra^j , (15) 

so that when the number of particles in the system Ej(nj) = X]i('^f) ~ is odd, the 
equivalent fermionic Hamiltonian satisfies periodic boundary conditions; whereas, if A^b 
is even, antiperiodic boundary conditions are required. 

The above mapping to noninteracting fermions allows one to realize that the 
presence of an additional periodic potential opens gaps at the edges of the reduced 
Brillouin zones. This implies that new insulating phases appear at fractional fillings 
Hi = i/£, with i = !,...,£— 1, in addition to the insulating phase at n = 1 which is 
present in the absence of the superlattice. 

In the following, we address the question of what happens during the evolution of 
initially prepared insulating states when they are allowed to expand in a superlattice. 
For simplicity, we restrict our analysis to the case i = 2. (The generalization to larger 
values of i is straightforward.) To study these systems, we follow the exact approach 
already described in detail in Refs. [121 ESI El • 

For i = 2, a. band gap A = 2 A opens at A; = 7r/2a. The dispersion relation for the 
two bands reads 

e±{k) = ±^At^ cos^ {ka) + A\ (16) 

where stands for the upper band and '— ' for the lower one. One can then calculate 
the group velocity at each momentum as 

^ de±{k) ^ 2t^sm{2ka) 
dk v^4t2 cos2(A;a) + 
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which means that the wave vectors at which the maximum group velocity occurs satisfy 



cos^ (kma) 



1 + 4t£ _ . 

A2 



:i8) 



For these values of k, the density of states attains its minimum values. Hence, under 
the appropriate initial conditions, we should expect the quasi-condensates to emerge at 
these values of k rather than aX k = 7r/2, as in the absence of the superlattice {A = 0). 




Figure 8. Evolution of (a) density and (b) momentum profiles of 101 HCB's, initially 
prepared in an insulating state with density one, on 1000 lattice sites. The superlattice 
parameters are £ = 2 and A = 2t. The times are r = (V), 50h/t (□), 200h/t (Q), 
and mh/t (A). 



Since for i = 2 the only insulating phases occur at full filling and at half filling, we 
start studying the evolution of an insulating state initially prepared with one particle 
per lattice site, like in the previous sections and in Refs. ^21 El- The expansion of such 
a state with 101 HCB's is shown in Fig.|H| While the evolution of the density [Fig. |H^a)] 
is very similar to that in the absence of the superlattice ^21 E] , the evolution of the 
momentum distribution function is completely different [Fig. IHl^b)]. No peaks appear in 
Uk, in contrast to the case analyzed in Refs. fT^ IT^ and in the previous sections. This 
is because in the superlattice the mean energy per particle (e = 0) for the fully filled 
insulating state lies in the band gap, and the states with the closest energy are the ones 
with momentum k = ±7r/2, which have i^^ = and the maximum density of states, i.e., 
exactly the opposite of the case without the superlattice. 

A scenario in the superlattice that is closer to the one in a Mott insulator with 
n = 1 (and finite U) in the absence of the superlattice, is the one in which the initial 
state is prepared with a mean density of 0.5. Such state has short-range (exponentially 
decaying) one-particle correlations like the ones seen in Fig. In addition, its mean 
energy per particle lies within the lowest band, where the density of states is finite. 
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Figure 9. Evolution of (a) the mean density per unit cell and (b) momentum profiles 
of 100 HCB's, initially prepared in an insulating state with a mean density of 0.5, on 
1000 lattice sites. The superlatticc parameters are t = 2 and A = 2t. The times are 
T = (V), bm/t (□), 2mh/t (O), and 400;i/t (A). 



In Fig. ini we show the evolution of the mean density per unit cell and momentum 
profiles during the expansion of a state prepared in a half-filled box with 100 HCB's 
and A = 2t. As can be seen in Fig. Efb), the initial momentum distribution is not fiat 
like the one in Fig. IHb). Its maximum ai k = signals the presence of short-range 
correlations like in the Mott insulator of Fig. |3J Remarkably, during the expansion of 
this state, sharp peaks appear in at ka = ±0.87 and ka = ±2.27, which are the 
momenta for which the group velocity is maximum and the density of states has a 
minimum, following Eq. (fTH|) . 

The peaks in signal the emergence of quasi-condensates of HCB's at finite 
momentum. This can be seen by studying the natural orbitals which can be 

considered to be effective single-particle states in interacting systems, and are defined 
as the eigenf unctions of the one-particle density matrix pij |41j . 

N 

J2p^^^(^W^) = ^r,{rmr), (19) 
i=i 

with occupations A^. In dilute higher dimensional gases, where only the lowest natural 
orbital (the highest occupied one) scales ~ A^^, the occupation of this orbital can be 
regarded as the BEC order parameter, i.e., the condensate 

In Fig. irUT a) we show the values of for the 200 highest occupied orbitals after 
the same expansion times as in Fig. IHl One can see that during the expansion the lowest 
natural orbital becomes "highly" populated with of the order of ^/Nb particles. This 
is because quasi-long range correlations develop in the system, as shown in the inset 
of Fig. ITUT a). The power- law decay of the one-particle correlations is ~ 1/ — Xj\, 




Figure 10. Evolution of (a) the natural orbital occupations and (b) the modulus 
of the lowest natural orbital wave function (averaged in the unit cell) for 100 HCB's 
initially prepared in an insulating state with mean density of 0.5 on 1000 lattice sites. 
The superlattice parameters are £ = 2 and A = 2t. The inset in (a) shows the 
initial exponential decay of one-particle correlations (averaged per unit cell), and its 
conversion to a power-law decay in the region where the quasi-condensate forms. The 
black line corresponds to 1/ \/\ Xi — Xj \. The times are t = (V), 50h/t (□), 2QQh/t 
(O), and 400h/t (A). 



i.e., the same form that appears in the absence of the superlattice f2l ^^"^ that 
has been proven to be universal in the ground state [T^ [T7j . The wave function of 
the lowest natural orbital during the expansion is depicted in Fig. ITOf b). One can see 
that it exhibits exactly the same features observed in Refs. [T^ IT^ when there was no 
additional periodic potential. After the Mott insulator melts, the shape of the lobes of 
the natural orbital stops changing and they just move with the maximum velocity in 
the lattice given by Eq. (fT7|) for k = km- 

One final remark on these emerging quasi-condensates is in order. In contrast to 
those emerging in a system without a superlattice potential which are mainly formed 
by particles with ka = ±7r/2, we find that in the superlattice the quasi-condensates are 
mainly formed by HCB's with the four momenta km- This is reflected by the four-peak 
structure of the Fourier transform of (p^ at momenta km, depicted in Fig. ^2 Hence, the 
four peaks that appear in the momentum distribution in Fig. IHl^b) reflect the formation 
of the quasi-condensates in Fig. [TUT b). which have a richer structure in /c-space than 
those that emerge in the absence of the superlattice. 

5. Conclusions 

In this work, we have presented non-trivial generalizations of a nonequilibrium 
phenomenon originally found in systems of hard-core bosons that expand out of strongly 
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ka 

Figure 11. Evolution of the Fourier transform of the lowest natural orbital wave- 
function of 100 HCB's, initially prepared in an insulating state with mean density of 
0.5 on 1000 lattice sites. The superlattice parameters are £ — 2 and A = 2t. The times 
are 50h/t (□), 200;i/t (Q), and 400^^ (A). 

correlated Mott-insulating states. First, we have shown that the emergence of coherent 
matter waves at finite momenta persists away from the hard-core hmit, i.e., to finite 
values of the on-site repulsion U which range down to the critical values, as long as 
a Mott-insulating state is attainable for the initial state. The accurate treatment of 
this complicated many-body problem has been made possible by the advent of the 
t-DMRG. By comparing to exact results in the hard-core case, we have shown that 
our t-DMRG calculations, although approximate, are very reliable. Two new features 
appear in the soft-core case: (i) Although the time evolution of the density profiles 
is indistinguishable from those of hard-core bosons at large values of U, the initial 
momentum distribution functions are markedly different. Nevertheless, after the Mott 
region melts, the momentum distribution functions of both systems become nearly 
identical, (ii) As the strength of the interaction is reduced, a shift of the momentum 
of the coherent matter wave to values smaller than 7r/2a is observed. The appearance 
of a power law in the spatial decay of the one-particle density matrix demonstrates 
explicitly the (quasi-)coherent nature of the resulting matter wave. Furthermore, we 
have shown that a coherent matter wave can be also obtained in the presence of 
superlattice potentials given an adequate selection of the initial insulating state. In 
the latter case, the emerging quasi-condensates have a more complicated structure in 
momentum space, with a number of dominating momenta whose locations depend on 
the superlattice used. 

The results of this work show that it is possible to engineer atom lasers with 
a high degree of control. The momentum of the coherent matter wave can first be 
regulated by setting the wavelength of the underlying optical lattice and further fine- 
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tuned by regulating the depth of the potentials in the lattice, i.e., the intensity of the 
corresponding laser beam. A further finite shift of the momentum can be achieved 
by superimposing another laser beam with a commensurate wavelength, giving rise 
to a superlattice, which leads to the emergence of a matter wave with two (or more) 
dominating momenta. 
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